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Instabilities of Spiral Shocks I: Onset of Wiggle Instability 
and its Mechanism 
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ABSTRACT 

We found that loosely wound spiral shocks in an isothermal gas disk caused by a 
non-axisymmetric potential are hydrodynamically unstable, if the shocks are strong 
enough. High resolution, global hydrodynamical simulations using three different nu- 
merical schemes, i.e. AUSM, CIP, and SPH, show similarly that trailing spiral shocks 
with the pitch angle of larger than ~ 10° wiggle, and clumps are developed in the 
shock-compressed layer. The numerical simulations also show clear wave crests that 
are associated with ripples of the spiral shocks. The spiral shocks tend to be more 
unstable in a rigidly rotating disk than in a fiat rotation. This instability could be 
an origin of the secondary structures of spiral arms, i.e. the spurs/fins, observed in 
spiral galaxies. In spite of this local instability, the global spiral morphology of the gas 
is maintained over many rotational periods. The Kclvin-Hclmholtz (K-H) instability 
in a shear layer behind the shock is a possible mechanism for the wiggle instability. 
The Richardson criterion for the K-H stability is expressed as a function of the Mach 
number, the pitch angle, and strength of the background spiral potential. The crite- 
rion suggests that spiral shocks with smaller pitch angles and smaller Mach numbers 
would be more stable, and this is consistent with the numerical results. 
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1 INTRODUCTION 

Spiral arms are the most prominent substructure in disc 
galaxies. The spiral density wave hypothesis (Lin & Shu 
1964) has been widely accepted as explaining the steady stel- 
lar spiral structure. The spiral density waves cause shocks in 
the interstellar medium, i.e. the galactic shocks. This phe- 
nomenon was well studied in the 60s and 70s. Compression 
of the interstellar medium behind the shocks is an essen- 
tial mechanism in the formation of massive stars that il- 
luminate the spiral structure. Assuming tightly wound spi- 
rals and steady flow, linear and weak nonlinear response of 
the gas in the spiral potential has been studied. (Fujimoto 
1968; Roberts 1969; Shu et al. 1972; Shu, Milione, & Roberts 
1973). 

Most of these initial studies adopted rather simple as- 
sumptions for the ISM, i.e. isothermal equation of state 
and single-phase, uniform and non-selfgravitating gas. Then, 
more realistic cases were studied. For example, Sawa (1977) 
investigated effects of the selfgravity of the gas on the 
asymptotic stational solution. Ishibashi & Yoshii (1984) 



* E-mail: wada.keiichi@nao.ac.jp 

f JSPS fellow, E-mail: koda@astro.caltech.edu 



studied an energy loss flow. Tomisaka (1987) took into ac- 
count the thermal heating and cooling processes for the 
cloud fluid. Lubow, Cowie, & Balbus (1986) studied two- 
component flow, i.e. isothermal gas and stars. Among these 
studies, the spiral potential was assumed to be tightly 
wound, and a periodic boundary condition was applied. 
These studies focused on the steady state of the gas in a 
spiral potential, because the observed spirals are thought to 
be a long-lived structure in galaxies. Woodward (1975), on 
the other hand, studied gas flow in a tightly wound spiral 
potential, using time-dependent numerical simulations un- 
der the same approximation adopted by Roberts (1969), and 
showed that shocks are formed within one or two transits of 
the gas through the spiral pattern. Johns & Nelson (1986) 
performed numerical simulations on a time-dependent, two- 
dimensional gaseous response in a spiral potential, and found 
that quasi-steady spiral shocks are formed. 

The next question concerns the stability of the galactic 
shock: Nelson & Matsuda (1977) showed that tightly wound 
spiral shocks are dynamically stable. Dwarkadas & Balbus 
(1996) discussed the linear stability of the spiral shocks as- 
suming a flat rotation curve and concluded that the flow 
is stable. Mishurov & Suchkov(1975) showed that the flow 
obtained by Roberts (1969) turns out to be stable behind 
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the shock front. On the other hand, it is unstable ahead 
of the shock. Recently, Kim & Ostriker (2002) performed 
two-dimensional, shearing box simulations of the magne- 
tized gas in a spiral potential in an effort to understand 
the origin of the secondary structure in spiral arms, namely 
'spurs' (Elmegreen 1980). They found that gaseous spurs are 
formed as a consequence of magneto-gravitational instabili- 
ties inside spiral arms (see also Balbus (1988) for the linear 
analysis) . Gravitational stability of the shocked gas in spiral 
arms has been well studied in the context of the formation 
of the giant molecular clouds (e.g. Elmegreen 1994). 

In most of the studies mentioned above, spiral shocks 
are assumed to be tightly wound, i.e. their pitch angles are 
smaller than 10°, and a flat rotation curve which is in gen- 
eral suitable for the galactic outer disk is also assumed. In 
this paper, we investigate stability of the gas flow in more 
general spiral potentials, with small (~ 5°) to large (~ 20°) 
pitch angles in various rotation curves. We perform two- 
dimensional, time-dependent, global hydrodynamical simu- 
lations of the gas in spiral and barred potentials using three 
high-resolution numerical schemes based on different con- 
cepts. We show that the spiral shocks generated in a non- 
self-gravttating thin disc may be dynamically unstable, in 
the sense that the shock front wiggles, and eventually knots 
are formed in the shock-compressed gas (section 2 and sec- 
tion 3). The physical origin of this instability in terms of 
the Kelvin-Helmholtz (K-H) instability in a shocked layer, 
as well as other implication from the numerical results will 
be discussed in section 4. A brief summary and conclusions 
are given in section 5. In Appendix, the Richardson criterion 
behind the spiral shock is approximately derived. 

Non-linear and long-term evolution of the instability 
and comparison with observations, such as spurs/fins in 
galaxies, and dynamics of the multi-phase, self-gravitating 
interstellar medium in the spiral potential will be discussed 
in subsequent papers. 



2 NUMERICAL SIMULATIONS 
2.1 Methods 

As we will see in section 3, hydrodynamical simulations show 
that the spiral shocked layer is unstable in some situations. 
In order to ensure that the instabilities are not caused by 
numerical artifacts, we apply three hydrodynamical schemes 
based on different physical and numerical concepts to the 
same models, i.e. an Euler mesh code (AUSM: Advection 
Upstream Splitting Method), mesh-based semi-Lagrangian 
code (CIP: Cubic-polynomial Interpolation), and particle- 
based Lagrangian code (SPH: Smoothed Particle Hydro- 
dynamics). We find that the instability commonly occurs 
among results with these three codes. All these schemes are 
standard methods in computational fluid dynamics and as- 
trophysical gas dynamics. We briefly describe the three nu- 
merical schemes below: 



comparison with Flux-difference splitting schemes (e.g. Roe 
splitting) and PPM (? ). The two-dimensional and three- 
dimensional versions of the scheme with a Poisson equation 
solver were applied to various problems of dynamics of the 
interstellar medium (Wada 2001; Wada & Norman 1999, 
2001, 2002; Wada, Spaans, & Kim 2001; Wada, Gerhardt, & 
Norman 2002). We summarize the essential points of AUSM 
below. 

The basic equations for the isothermal two-dimensional 
flow are written as 
dU dF dG „ 

where U T = (p, pu, pv), F T = [pu, pu 2 + p, puv], G T = 
[pv,pvu,pv + p]. 

Flux F and G consist of two physically distinct parts, 
namely advection and pressure terms: 





v + 




(2) 



In AUSM, these two terms are separately split at a cell sur- 
face. The van Leer-type flux splitting process is applied sep- 
arately to these two terms (Liou 1996). 

We achieve third-order spatial accuracy with the Mono- 
tone Upstream-Centered Schemes for Conservation Laws 
(MUSCL) approach (? ). To satisfy the TVD (Total Vari- 
ation Diminishing) condition using MUSCL, we introduce 
a limiting function. The code was checked for various test 
problems, such as one-dimensional and two-dimensional 
shock tube problems, collision of strong shocks, reflection 
of shocks at a wall, and propagation of 2-D blast waves. See 
details in Wada & Norman (2001). We use a Cartesian grid 
with 128 2 - 2048 2 zones. 



2.1.2 CIP 

The Cubic-polynomial Interpolation or Cubic Interpolated 
Propagation (CIP) method, developed by Yabe et al. 
iYabe fc Aoki 199ll) . is a kind of semi-Lagrangian method, 
in which the hyperbolic equation df /dt + (u ■ A)/ = g, 
where e.g. / = (p,u) and g = (— pA ■ u,—Ap/p), is split 
into advection and non-advection phases. The former is ad- 
vanced by using a profile interpolated with a cubic polyno- 
mial, and the profile is determined from the time evolution 
of the quantity / and the spatial derivative A/. The CIP 
method is a stable and less diffusive scheme for compressible 
and incompressible flows with sharp boundaries, and it has 
been applied to various problems in the field of computa- 
tional fluid dynamics and astrophysics (e.g., Kudoh, Mat- 
sumoto & Shibata 1998 for formation of astrophysical jets) . 
One remarkable feature of this code is that discontinuity is 
well resolved with relatively small numbers of grid points. 
The code used here has first-order accuracy in space and 
time. The number of grid zones used here is the same as 
that which runs with AUSM. 



2.1.1 AUSM 

AUSM was developed by Liou & Steffen (1993), and has 
been widely used for aerodynamical simulations. AUSM is 
remarkably simple, but it is accurate enough based on a 



2.1.3 SPH 

We use the Smoothed Particle Hydrodynamics (SPH) code 
presented by Koda & Wada (2002), which basically uses 
Benz's formulation (Benz 1990), using the spline kernel by 
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Monaghan & Lattanzio (1985) with the modification for its 
gradient (Thomas & Couchman 1992). The correction term 
for viscosity is taken into account to avoid large entropy gen- 
eration in pure shear flows (Balsara 1995). The SPH smooth- 
ing length h varies in space and time, keeping the number 
of particles within the radius of 2h at an almost constant 
of 32 according to the method of Hernquist & Katz (1989). 
The leapfrog integrator is adopted to update positions and 
velocities. 

In order to represent an initial uniform-density gas disc, 
we randomly distribute 10 s gas particles in a 2-D disc with 
a radius of 3 kpc. The gas particles follow pure circular ro- 
tation that balances the centrifugal force. 

2.2 Basic equations and Boundary conditions 

We solved the following equations numerically in two- 
dimensions: 
3E 



^+V.(E 9 .) 
dv , „. Vp 



o, 



(3) 
(4) 



where, E 9 ,p, v are surface density, pressure, and velocity 
of the gas. 3> ext is a fix external potential. We assume the 
isothermal equation of state, and gas temperature is 10 4 K, 
unless otherwise mentioned. 

In AUSM and CIP simulations, the boundary condi- 
tions are set as follows. In ghost zones at boundaries, the 
physical values remain at their initial values during the cal- 
culations. From test runs, we found that these boundary 
conditions are much better than other boundary conditions, 
such as 'outflow' boundaries, which cause strong unphysi- 
cal reflection of waves at the boundaries. Even though this 
boundary condition does not seem to affect the gas dynam- 
ics in the computational box, we put the boundaries far 
enough from the region where we study the stability of the 
spiral shocks to avoid unexpected numerical artifacts. The 
size of the computational box is typically 8 kpc x 8 kpc 
(model A, B, and D, see below) or 32 kpc x 32 kpc (model 
C). In SPH, we use a free boundary. The simulations are 
performed in a rotating frame of a fixed pattern speed of 
the non-axisymmetric potential. We assume initially uni- 
form surface density (7M@ pc -2 ), 

2.3 Potential models 

We solve two-dimensional gas dynamics in a time- 
independent, non-axisymmetric external potential. Four dif- 
ferent potential models are used; three spiral potentials 
(model A, B, and C) and one bar potential (model D). One 
of the spiral models (model C) is assumed a large scale disk 
(initial radius is 16 kpc), and others and the bar model are 
for inner galactic disks (initial radii are 4 kpc). Rotation 
curves derived from the axisymmetric potentials of model 
A, B, and D, are shown in Fig. Q and model C in Fig. 
The four models are defined as below: 



model A 
model B 
model C 
model D 



$ cxt = MR) + MR) + MR, 
$ cxt = MR) + MR) + MR, 
$ cxt ee MR) + MR) + MR, 



(5) 
(6) 

(7) 
(8) 



Table 1. Parameters for the four potential models. Units of the 
core radii are kpc, and km s — 1 for the velocity v a , v\,, and v c . 



model 


a 


b c 


Va 


Vb v c 


A 


1.4 




220 




B 


0.2 


1.7 


150 


103 


C 


1.3 


6.8 


220 


250 


D 


1.4 


0.01 


220 


311 



where the axisymmetric components $0, $2, and $4 are de- 
fined as 



MR) = av 2 a (27/4) 1/2 (R 2 + a 2 )- 1/2 , 
MR) = K 2 (27/4) 1/2 (J? 2 + 6 2 )- 1/2 , 
MR) = cu c 2 (27/4) 1/2 (i? 2 + c 2 )- 1/2 , 



(9) 
(10) 

(11) 



and the constants a, b, c, v a ,Vb and v c in the four models are 
summarized in Table 1. The non-axisymmetric components 
(i.e. spiral and bar potentials), $1 and $3 are given by 



MR,< 



MR,' 



£0 



aR 2 § 







£<>- 



(R 2 + a 2 ) 3 /2 

bR 2 $2 



cos[20 + 2 cot i ■ ln(R/R(()\2) 



cos[2<£ + 2 cot i ■ ln(i?/#o(P) 



(R 2 + 6 2 ) 3 / 2 

where i is the pitch angle of the spiral potential (see Fig. ll2t . 
and i?o is an arbitary constant (Rq — 0.9 kpc is used for all 
models). A constant eo represents maximum strength of the 
non-axisymmetric potential. For the bar potential model D, 
i = 7r/2 is assumed with $1. 

Model D is chosen in order to see stability of spiral 
shocks formed by a different mechanism, i.e. resonance- 
driven spiral shocks. We add a point mass at the galactic 
centre ($4), which resembles, for example, a supermassive 
black hole. In model D, the 'black hole' mass is assumed 
to be 1% of the total dynamical mass. This causes a 'nu- 
clear' inner Lindblad resonance (nILR), and it is known that 
trailing spiral shocks are formed around the nILR (Fukuda, 
Wada, & Habe 1998). This characteristic resonant-driven 
structure also appears around the outer ILR (Athanassoula 
1992; Piner et al. 1995). Using two-dimensional SPH sim- 
ulations, Fukuda, Habe, & Wada (2000) revealed that the 
trailing spirals can be gravitationally unstable, and claimed 
that this process is important in fuelling the gas to the galac- 
tic centre. Therefore, stability of such trailing spirals in the 
galactic nuclear region in cases in which the gas is not self- 
gravitating is worth studying. We can also study the gener- 
ality of the stability/instability of the spiral shocks. 



3 NUMERICAL RESULTS 

In Fig. EH we show typical numerical results obtained by 
the three different numerical schemes, i.e. SPH, AUSM, and 
CIP, explained in section 2.1. 'Wiggle' instabilities of the spi- 
ral shock, namely ripples of the shock fronts and less dense 
fin-like structure are clearly seen in all the results. The global 
morphology of the spirals in all results is quite similar. Com- 
paring Fig.|3Ja) with^c), one may find that the instability 
is less prominent in the result by CIP. This is because the 
CIP code used here is less accurate (i.e. first-order scheme) 
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250 




0.0 0.5 -.0 1.5 2.0 

Radius (kpc) 

Figure 1. Rotation curves adopted in the numerical simulations. 
Solid line (model A): nearly rigid rotation models, dashed line: 
nearly flat rotation curve (model B), and dotted line (model D): 
the bar model with a central massive black hole. 



300 




Radius (kpc) 

Figure 2. The rotation curve of model C. Solid line: $o + 3>2, 
dotted line: <3?o, an d dashed line $2- 



than AUSM (spatial third-order scheme), and therefore the 
shock is captured less sharply. Although the SPH result is 
noisier and seems to already be in a non-linear phase espe- 
cially in the inner region, comparing to the other two results 
obtained by the mesh-based codes, wavelengths of the insta- 
bility are similar between the models. The fin-like structures 
seen in the inter-arm regions are density waves originating in 
the wiggle instability of the shock front. This phenomenon 
will be discussed in detail in a subsequent paper (see also 
section 5). 

For further discussion, we use results with the AUSM 
scheme. The four panels of Fig. 2] show how the numerical 
results are affected by the spatial resolution. Although global 
morphology of the spirals is similar in all results, the wiggle 
instability of the shock fronts is not clear in the results with 
128 2 and 256 2 cells. We need at least 512 2 grid points to 
resolve the instability in this case. Fig. [I] (d) shows that 
characteristic scale of the ripples is roughly the width of 




x(kpc) 

Figure 3. Response of the gas disc to a spiral potential (model A 
with i = 10°, £ = 0.1, n p = 26 km s" 1 kpc" 1 ) at t = 24 Myr. (a) 
Density distribution given by AUSM with 1024 2 cells. Gray-scale 
represents log-scaled surface density with a unit of Mq pc — 2 . (b) 
SPH with 10 s particles (roughly there are 4 X 10 4 particles in this 
box). Gray-scale represents amplitude of the spiral potential, and 
the white line is a bottom of the potential, (c) Same as (a), but 
by CIP. 



the shock-compressed layer, which is about 100-200 pc. We 
can also see that spur-like structures are developed into the 
inter-arm regions, which are associated with the ripples (see 
also Fig. 

From the azimuthal density profile (Fig. , it is clear 
that reveals sharp shocks where density jump is ~ 100, 
and large density fluctuation in the shock compressed layer, 
which is originated from the clumps seen in Fig. 0] are 
formed. One may also notice density fluctuation in the inter- 
arm regions. This can be seen as many 'crests' in Fig.0](d). 

Figure |S| shows time evolution of the spiral shocks in 
the same model of Fig. |3](d). The wiggle instability begins 
in the inner region first, then the outer shocks becomes un- 
stable. At t — 28 Myr (d), the instability is in a non-linear 
stage, and 'spurs' are developed in the inter-arm regions 
towards the other spiral shocks. Typical time-scale of the 
linear growth of the instability is ~ 10 Myr at R ~ 500 pc 
and ~ 30 Myr at R ~ 2 kpc. This is about a half of the or- 
bital time-scale. The spurs are nearly regulary spaced, and 
their intervals seem to be larger in the non-liear phase than 
in the early phase. 

The two panels in Fig.Qdemonstrate that the stability 
of the shock fronts is sensitive to the pitch angle of the spiral 
potential (model A). It shows that the shocks are stable for 




Figure 4. Effect of changing a spatial resolution on the instabili- 
ties. The potential model is the same as that in Fig- EI ( a ) Density 
distribution of a central 1/4 region of the computational box with 
128 2 cells at t = 24 Myr. (b), (c) and (d): Same as (a), but with 
256 2 , 512 2 , and 2048 2 , respectively. The parameters used here are 
the same as those for models in Fig. 3. 
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Figure 5. Azimuthal density profile at R = 1.5 kpc of Fig.0J(d). 
The unif of density is Mq pc -2 . 



i = 5°, whereas they are unstable for more loosely wound 
spirals, i.e. i = 20°. 

We found that if the rotation curve is fiat (model B), 
the spiral shocks are stable (Fig-EJ, provided that the spiral 
potential is weak (but strong enough to produce the shocks). 
Note that even if the unperturbed rotation curve is fiat, 
the perturbed potential can cause a radial gradient in the 
circular velocity on a local scale. This local velocity gradient 
is essential for the onset of the instability as discussed in 
section 4 and Appendix. 

We found that the wiggle instability of the shocked layer 
appears not only in the central kpc of galaxies, but also in a 
larger size disk. Fig.|5|is a snapshot of model C, in which a 
disk with four times larger size than that in model A and B 
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Figure 6. Evolution of the instability in the same model shown 
in Fig.[I](d). Four snapshots are at (a) t=7.5 Myr, (b) 12.5 Myr, 
(c) 17.5 Myr, and (d) 28 Myr. 




xftoc) x(kpc) 

Figure 7. Response of the gas disc to a spiral potential (model 
A with i = 5° and i = 20°). e = 0.1, C p = 26 km s" 1 kpc~ x 
and T g = 10 5 K at 50 Myr. 1024 2 cells are used in both models. 




Figure 8. Same as Fig.UI but for the case that the rotation curve 
is nearly flat (model B) with eo = 0.03. 



6 Keiichi Wada and Jin Koda 




Figure 9. Density distribution in a spiral potential (model C) 
with i = 15° at t = 10 8 yr. 2048 2 cells are used. 



is used. The instability and spur-structure in the inter-arm 
regions are similar between the model A and C, suggesting 
that the wiggle instability can occur on any scales. 

Finally, we examine whether the instability seen in the 
previous section is a phenomenon peculiar to the spiral po- 
tential. Fig. 1101 is the density and velocity field near the 
nuclear Lindblad resonance caused by a central mass con- 
centration, such as a supermassive black hole, in a weak bar 
potential (model D). In this case, off-set straight shocks are 
formed outside the nuclear spiral shocks. We found that the 
shocks are unstable, and they wiggle as seen in the snap- 
shot at t = 70 Myr. One may also notice the wavelets in the 
post-shock regions. 



4 DISCUSSION 

The results shown in section 3 strongly suggest that the 'wig- 
gle instability' of the shocked layers followed by formation 
of clumps and spurs is not caused by numerical artifacts, 
but it has a physical origin. We suspect that they originate 
in the Kelvin-Helmholtz (K-H) instability. As the velocity 
fields in Fig. llOl and llll show. there is a strong velocity shear 
behind the shock. In standing spiral shocks caused in a rotat- 
ing disk, like galactic shocks or bar-driven shocks, it is most 
likely that velocity gradient is generated behind the shocks. 
In general, the K-H instability can develop in a contact sur- 
face of two fluid layers with different tangential velocities, 
even if the velocity difference is small, provided that there 
is no surface tension. However, in the situation considered 
here, buoyancy force can suppress the instability, because 
there is also a density gradient behind the shock under the 
spiral (or bar) potentials. The Coriolis force could also affect 
development of the K-H instability (Chandrasekhar 1961). 
However, for the sake of simplicity, we ignore this kind of 
effect in the following discussion. 

A condition of the K-H stability in a sheared layer with 
a density gradient is characterized using a non-dimensional 
number, i.e. the Richardson number, J, defined by 



S 




-1.0 



Figure 10. Spiral shocks generated around the nuclear Lindblad 
resonance. ( model D with eo = 0.1, C p = 30 km s — 1 kpc - 1 ) at 
t = 60 (upper panel) and 70 Myr (lower panel). Gray-scale is the 
log-scaled density, and vectors represent a velocity field. 



where g is gravitational acceleration, which is normal to the 
shock front (^-direction), and u is shear velocity. J mea- 
sures the ratio of the buoyancy force to the inertia, and a 
necessary condition for the K-H stable is that J > 1/4 (e.g., 
Chandrasekhar 1961) . Here we give a rough estimate of J be- 
hind the shock. Figure [T2*l shows flows passing a spiral shock 
schematically. The K-H instability could be originated in the 
velocity shear (Am) between position B and C, i.e. v" (> v') 
and v. Equation 1141 is then J ~ 2gd/Au 2 for the shocked 
layer with a width of 2d. Then using 



9 ■ 



Rdn 



ev. 



4> 



Rs'mi 



(15) 



where $ na is a non-axisymmetric component of the potential 
(e.g. $i or $3 in section 2.3), e represents a fraction of the 
non-axisymmetric part to the axisymmetric potential, and 
t)0 is circular velocity, and 



J 



g dpjdz 
p (du/dz)' 



(14) 



Au « (v" - v) cos ittv+fl- j^j 



COS I, 



(16) 
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Figure 11. Velocity distribution in the early phase of the evo- 
lution for model A with i = 10°, eo = 0.1. Gray-scale represents 
amplitude of the velocity. 
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Figure 12. Spiral shock and streamlines, for which we assumed 
that the curvature can be negligible, and therefore the spiral 
shocks are modeled as oblique shocks. For AR -C Ft, the spi- 
ral shock can be approximated as an oblique shock. A streamline 
is changed its direction at point A with velocity v' , and it is ac- 
celated toward the other shock. At point B, its velocity is v" . A 
velocity gradient is expected between points B and C. A char- 
acteristic scale of the shocked layer is 2d. (£,??) is a orthogonal 
coordinate along the spiral with the pitch angle, i. 



where M is the pre-shock Mach number, the Richardson 
number can be expressed as a function of the Mach number 
and the pitch angle: 

2ed 



J 



i?sini(l-sin 2 i)(l-l/M 2 ) 2 ' 



(17) 



This suggests that for a given M, a larger J is obtained for 
smaller i, and for a larger M, J is smaller. For example, for 
e = 0.1 , M = 10, and d/R = 0.1, J ~ 0.12 and 0.39 for 
i — 10° and i = 3°, respectively. Although one cannot pre- 



dict whether the flow is unstable only from the Richardson 
criterion, eq. 1171 suggests that tightly wound spiral shocks 
may be relatively stable, comparing to open spirals. This 
tendency is in fact confirmed by our numerical simulations 
(e.g. Fig. [7J . More detail treatment of the flow and deriva- 
tion of the Richardson number behind shocks is given in 
Appendix, where Fig. H3l shows how J depends on the Mach 
number and the pitch angle, which is qualitatively the same 
sense of eq. 1171 1. 

As the ripples develop in the shock compressed layers, 
many spur-like structures, or "crests", which are associated 
with the clumps, are formed behind the shock, and they 
extend to the inter- arm regions with a large pitch angle. 
These structures are morphologically similar to the dust- 
lanes observed in spiral galaxies. A typical example of this 
is the central regions of M51 (Scoville et al. 2001, see also 
Elmegreen 1980). Kim & Ostriker (2002) argued that those 
spur-like structures are formed as a consequence of magneto- 
Jeans instability using two-dimensional, local MHD simu- 
lations. The wiggle instability we found in the numerical 
simulations would be an another mechanism to form the 
spurs, which originates in the K-H instability of the shock- 
compressed layer. The mechanism of forming spurs is follow- 
ing. The clumps formed by the instability have an internal 
velocity/angular momentum gradient that causes a phase 
shift in the orbital motion, and results in deformation of the 
clumps. Hence, a part of each clump with relatively small 
angular momenta falls along the spiral shocks, whereas the 
gases with larger angular momenta try to maintain nearly 
circular rotation. As a result they move away from the shock 
toward the other spiral shock. Because of the internal gra- 
dient of the angular momentum of a clump formed behind 
a shock, the clump gas eventually extends to the inter-arm 
region as it rotates in the galactic potential. One should 
note that even in a highly nonlinear phase, the global galac- 
tic "shocks" are still present, although they are no longer 
continuous shocks. In subsequent papers, we will discuss 
the non-linear development of the ripple instability of spiral 
shocks and the formation mechanism of the spurs. Detail 
comparison with observations would be also necessary. 

In the present simulations, we ignored self-gravity of the 
ISM. Self-gravity is essential for dynamics and structure of 
the ISM in galaxies, especially for forming high-density re- 
gions and star formation. This is more prominent under the 
non-axisymmetric external potential (Wada & Koda 2001). 
As shown in Fig. |S] spiral shocks produce compressed lay- 
ers where the density is a factor 100 larger than that in 
the pre-shock region. This is because we assumed isother- 
mal equation of state for the ISM, and the Mach number is 
~ 10. Although this density contrast is comparable to a ra- 
tio between a typical density of giant molecular clouds and 
the average ISM density in a galactic disk, inclusion of self- 
gravity with realistic cooling processes is necessary to discuss 
hydrodynamical instabilities in real galaxies. Moreover, it is 
of interest to study how the gravitational instability of the 
shocked layer (Balbus & Cowie 1985 and Balbus 1988) re- 
lates with the wiggle instability. Even if the post-shock layer 
is gravitationally stable, self-gravity would trigger collapse 
of clumps formed by the wiggle instability, and it could lead 
star formation. 

Another interesting phenomenon related with the wig- 
gle instability is interstellar turbulence. In a non-linear 
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phase, clumps and spurs interact each other, and also collide 
with other spiral shocks. As a result, irregular motion of the 
gas is generated, and it could develop turbulent motion of 
the ISM. 

One might think that our findings seem inconsis- 
tent with previous studies on spiral shocks. For example, 
Dwarkadas & Balbus (1996) claimed that the spiral shocks 
are K-H stable, using a linearized perturbation analysis of 
the post-shock flow. In their analysis, they numerically in- 
tegrate a set of linearized perturbed equations assuming a 
tightly wound spiral (Roberts 1969). They found that the 
amplitude of the perturbation increases initially for about 
one-third of a rotation period, and then it decreases by a 
factor of 20-50 from its maximum value toward a steady 
state within one rotation period. It should be noted, how- 
ever, that they assumed a tightly wound spiral in a flat 
rotation curve (k = 2fi), which is a preferable condition 
for the K-H stable as shown in the numerical results and 
analysis in Appendix. The other point is that one-third of 
a rotation period would be long enough to form the rip- 
ples, clumps and spurs. These sub-structures found in our 
numerical simulations are not steady, but rather, temporal 
from a Lagrangian point of view. However, successive for- 
mation and destruction of these substructures could keep 
global, long-lived non-axisymmetric, non-uniform features 
in the interstellar medium in spiral galaxies. 

The instability that we found is not apt to appear in 
tightly wound spirals with a flat rotation curve. Therefore 
it is not surprising that the instability has not been found 
in the many studies in the 60s and 70s on spiral shocks. 
Yet there have been some pioneering studies on the time- 
dependent, two-dimensional flow in a spiral potential with 
various pitch angles. Johns & Nelson (1986) performed such 
simulations using two distinct codes, i.e. a Eulerian grid code 
and SPH. Unfortunately their spatial resolution was rather 
limited (63 2 x 2 polar grid), and as we showed in section 3, 
this is not fine enough to resolve the instability. However, in- 
terestingly, they noticed secondary structures in their spiral 
arms. The spiral arms in their simulations are not smooth, 
and local density peaks are seen in the contour maps. They 
mentioned a possibility of a physical origin of the modulation 
of the main two-armed response rather than some stochas- 
tic or numerical effect, however they did not confirm it. The 
secondary structure that they found might come from the 
instability of spiral shocks discussed in this paper. 



5 CONCLUSIONS 

Using high-resolution numerical simulations, we have found 
that spiral shocks caused by a galactic spiral potential and 
by bar-driven resonances can be unstable, in the sense that 
shocked layers are rippled, and in a non-linear phase they 
fragment into clumps followed by forming spurs. The insta- 
bility is most likely caused by the Kelvin-Helmholtz insta- 
bility in the post-shock layers where density gradient and 
velocity shear are present under the influence of the grav- 
ity of the spiral potential. This instability appears in spi- 
ral shocks on small (R < kpc) to large (~ 10 kpc) scales, 
also suggesting that local velocity shear behind the shocks 
drive it. Spiral shocks tend to be more unstable, if the their 
pitch angle is larger (> 10°), and for a higher Mach number 



(M > 3). The numerical results are not inconsistent with 
the discussion using the Richardson criterion (section 4 and 
Appendix). However, one should be careful that the Richard- 
son criterion itself is just a necessary condition for the K-H 
stability, and even if J < 1/4 the flow could be stable. In 
this sense, the K-H instability is one possible mechanism of 
the wiggle instability. We would need careful linear analyses 
for perturbations in loosely wound spiral shocks for further 
discussion. 
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APPENDIX: RICHARDSON CRITERION 
BEHIND A SPIRAL SHOCK 

In this Appendix, we stand on a hypothesis, i.e. the wig- 
gle instability found in the numerical experiments is caused 
by the Kelvin-Helmholtz instability. We give the Richardson 
number behind the spiral shocks as a function of parameters, 
such as the Mach number and the pitch angle of the spiral 
potential. This is useful to evaluate whether the sheared flow 
can be Kelvin-Helmholtz stable. One should note, however 
that we simplify the kinematics of the gas in spiral poten- 
tial, and the Richardson criterion here is not derived from 
a dispersion relation in a sheared layer caused by open spi- 
ral shocks, which may be ultimately given in a future work 
(see Dwarkadas & Balbus (1996) for linear analysis of tightly 
wound spiral shocks in a flat rotation curve). Yet a simple 
'toy model' below would be useful to understand physics 
behind the numerical results. 

We simply assume almost straight streamlines behind 
the shock and AR/R <C 1 (Tig. 1121 as we did in section 4. We 
also assume that the flow is isothermal, and the asymptotic 
solutions of the flow in a spiral potential is used to have the 
velocity and density gradient in the post-shock layer. First 
we consider only velocity shear which is parallel to the shock 
as a source of the K-H instability. Effect of the perpendicular 
velocity, v$, will be considered later. 

Under these assumptions, the Richardson number [eq. 
114t j can be expressed as 

J ~ god/ul, (18) 

with the velocity and density differences in a shocked layer, 
i.e. Mq = ( v c - «b) 2 /4 and a = (pc — Pb)/(pc + Pb) (B and 
C denote positions shown in Fig. 1121 , and 2d is width of the 
shock layer. 

Width of the shock layer, 2d, is related to the Mach 
number M (i.e. v^/cs sini, where v,/, is the pre-shock velocity 
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on a rotating frame with a pattern speed Q p ) and to the 
pitch angle i as 



2d 



AR 



tan i/M 



AR 



(19) 



sin i M 2 cos i 

As shown in Fig. 1121 a flow passing the oblique shock has 
velocity v' at position A, and is accelerated toward the 
other spiral shock. As a result, a velocity gradient is gen- 
erated behind the shock, for example between positions 
B and C (i.e. v^q and v'^ . Subscript "0" denotes unper- 
turbed variables). This velocity gradient, which is normal 
to the spiral shock, may be a source of the K-H instabil- 
ity. The post-shock velocity transverse to the shock, v v o, is 
v v o(R) = v$(R) sin i/M 2 , and the post-shock velocity par- 
allel to the shock is v^o{R) = v$(R) cosi, and the post- 
shock velocity at position A is v'(R) = \/ v2 q + 1,2 o = 

Vj,(R)*J cos 2 i + sin 2 % /M 4 . The post-shock velocity at po- 
sition B is 



''50 



(R + 2d) = v iQ (R + AR) + Avi, 



where A£ 
Av{. = 



= AR/(R sin i), and 
dv ( (R + AR) 



(20) 

(21) 
(22) 
(23) 

v v0 20 c M 2 cos i 

Here, we use dv^/d^ — dv^/dr) ■ Ar]/A£, with Arj/A^ = 
tan i/M 2 , and the linear solution of flow in a spiral potential 
(e.g. Spitzer 1978 [equation (13-24)]): 



(R + AR) 



V„l K 



-AC 

V V 1 K 2 

v v o 2Q, C 
AR 



tani AR 
M 2 i?sin? 



drj 



where fi, 
given by 



drj 



V V 1 K 



= v<j>/R + Op, 



(24) 
(25) 

v^sini/M 2 , and w^i/w^o is 



c 3 

V 2 



+ K* 



-v v a 



d 2 $ s 
R?dr] 2 



(26) 



with the epicyclic frequency k and a spiral potential $ s 
Using equation 1261 . equation 1231 is rewritten as 



Aw e 



AR 



y r 



) 2 (l 



+ K 2 



2tt c M 2 cos i 



R^dri 1 2n c M' 2 cosi 



[-2(0 + 1)-! (l - M 2 /sin 2 i) + l] 



(27) 



(28) 



where we assume a rotation curve as v$ = vo(R/ Ro) a , there- 
fore k 2 = 2(a + i)(vtf,/ R) 2 , and n^o = v^s'mi/M 2 is used. 
Using <3? s (r;) = — e$o(-R) sm(2-7rJ7/7rsini), in the vicinity of a 
spiral shock, 

#*. . ^ • (29) 



£?77 2 



4e$o . / 2r; 
— =— sin — ; — — 

sin i V sin % 



2A V 

sin i ' 



where Arj = Ai?/ (RM 2 cos i). For a slow pattern speed, i.e. 
v$ ~ i?fi c , from equations 1281 and 1291 then 

Au 5 



p<5 



(30) 



where 5 = AR/ R, and 

4e/(Af 4 cos 2 i • sin 3 i) 



p(M, i, a, e) 



I - 2(a+ l)- 1 (l - Af 2 /sin 2 i)] 



(31) 



We finally have the velocity variation «o m equation 1181 
using p an 5: 

[(v' ! ! (R + 2d)-v 60 (R))/2] 2 
4 



2 
''-() 



— |Ai?cosi + A« 5 

ail 



(32) 
(33) 



-p [ac5 cos i + /z<5 2 ] 



(34) 



where we use v^o(R + AR) — V(q(R) — [v$(R + AR) — 
v<t>(R)] cosi . 

The density gradient behind the shock, a, can be es- 
timated using the mass conservation, namely, poM 2 v'(R + 
AR) = p"v"(R + 2d), which is equivalent to 

p M 2 v^(R + Ai?)(sin 2 i/M 4 + cos 2 i) = 

a a I -. 

P «£0 I 1 + 



9 \ I/ 2 

tan i 
M 4 



(35) 



where v v /v^ = tan i/M 2 is used. Therefore, the dimension- 
less density change in the layer, a, is 



PoM 2 



p M 2 + p" 



1 — cos i ■ v$/ v'^q 

1 + cosi ■ V4,/v'^ 

a5 cos i + p,8 2 

2 cos i + fj.5 2 + aS cos i 



(36) 



(37) 



The gravitational acceleration, transverse to the shock 
front, g, is given by 



9 



<9<E> S 
Rdr] 



R 



2s 



25 



M 2 > 



(38) 



(39) 



Finally, the Richardson number, J, for the flow behind 
a spiral shock is given combining equations (1191 , 1341 , 1371 , 
and 1391 as 



J(i,M,e,a; A) 



4^ cos f j 2 ^) 

sin % \ cos 2 I 



M 2 cos i(a cos i + /xA sin i) (2 cos i + /xA 2 sin 2 i + aA cos i sin i) 

where A is a wavelength normalized by the radius, A = 
X/R — 8/ sini, and p, is given by equation l.'ilt . 

In Fig. 1131 the Richardson number given by equation 
1401 is plotted as a function of the Mach number, M, and 
the pitch angle, i. It shows that J monotonically decreases 
with M, therefore weaker shock are expected to be more 
stable. For a given M, a larger J is obtained for smaller i, if 
i < 20° . This suggests that tightly wound spiral shocks may 
be relatively stable. 

Equation 1401 is approximated as J « 4e/(iaM 2 ), for 
i < 1, A < 1 and a > 0. This suggests that spiral shocks in 
the region of a flat rotation (i.e. a — 0) curve are expected to 
be stable for short wavelengths in a galactic outer region. In 
fact, one of the numerical results (Fig.|5J shows that shocks 
in a flat rotation curve is stable. One should note, however, 
that even if the unperturbed circular velocity is constant, a 
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Figure 13. Richardson number, J, behind a spiral shock as a 
function of Mach number, M, and the pitch angle of the spiral, i, 
for a = 1 (i.e. rigid rotation ), A = 0.1, and strength of the spiral 
potential, e = 0.1. 




5 10 15 20 25 30 35 40 

i 



Figure 14. Curves for J = 1/4 are plotted as a function of the 
pitch angle i (degree) and Mach number M, below which the 
spiral shocks could be Kelvin-Helmholtz stable (i.e. J > 1/4, see 
Fig.EJand the text). Thick curves are for rigid rotation (a = 1), 
and thin curves are for Keplerian (a = —1/2) with e = 0.1 (solid 
curves) and 0.05 (dashed curves). A = 0.01 is assumed. 



local gradient in circular velocity caused by a spiral potential 
could destabilize the spiral shocks on a local scale. 

In Fig. 1141 the critical Richardson numbers are plotted 
as a function of the Mach numbers and the pitch angles for 
two different rotation curves. One should note again that 
J > 1/4 does not necessarily mean that the sheared layer is 
K-H stable. This is only a necessary condition for stability 
as mentioned in section 4. 

The curves in Fig. 1141 have minima, M m i n ~ 3 at pitch 
angles i m - m ~ 30. The plot shows that the J = 1/4 curves de- 
crease with the pitch angle, i, and they increase for i > i m i u . 
For i < imin, a flow is less affected by an oblique shock with 
smaller pitch angles, therefore the velocity gradient behind 
the shock is smaller, and as a result the shocked layer is 
expected to be more stable. On the other hand, if the pitch 
angle is larger than i m i n , the width of the sheared layer, 2d, 
becomes larger, and again the layer becomes stable. We can 
also know that weak shocks with M < M m i n are expected to 
be stable for all pitch angles. M m i n is larger in the Keplerian 
rotation curve than for the rigid rotation. This means that 
spiral shocks generated around a point mass (e.g. an accre- 
tion disk in a close-binary) tend to be more stable than those 
generated in a potential with a constant mass density, for 
the same pitch angle and the Mach number. 

For a given rotation curve, shocks may be more stable 
in stronger spiral potentials. This is because the K-H insta- 
bility can be suppressed by buoyancy force, and stronger spi- 
ral potentials cause larger gravitational acceleration nearly 
perpendicular to the shock front [see equation (1181 1. One 
should note however that Mach number of the flow towards 
the potential minimum is also affected by the strength of the 
potential, therefore it is not so straightforward to evaluate 
the effect of the spiral potential on the wiggle instability. If 
the spiral potential is too weak (e.g. amplitude of the spiral 
component is less than a few % of the axisymmetric one), 
shocks do not appear in the flow (e.g. Shu et al. 1973). 



In above discussion, we neglect the effect of expansion 
(perpendicular to the shock) velocity, v^. In fact, the flow 
behind the shock has the expansion velocity, which is ap- 
proximately 

Avr, a (t& - ^o)At7/A£ = 2u tani/M 2 = ( 4 1) 

where Mo = v^/c s . If Av^ > u^, we expect that the shear 
layer is stable for the K-H instability due to the expansion 
flow. From eq. (1411 . it is suggested that the expansion veloc- 
ity becomes important for spirals with a small pitch angle 
and/or small Mach number (see also Dwarkadas & Balbus 
1996). For example, this can be happen when Mo < 2.5, 
3.4, and 4.8, for the spirals with i = 20, 10, and 5°, respec- 
tively. In other words, for strong shocks with a large pitch 
angle, the post shock flow is nearly parallel to the shock, 
and then the expansion velocity on the K-H instability can 
be negligible. This is consistent with our numerical results. 
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